Collective states in a ring network of theta neurons

We consider a ring network of theta neurons with non-local homogeneous coupling. We analyse the corresponding continuum evolution equation, analytically describing all possible steady states and their stability. By considering a number of different parameter sets, we determine the typical bifurcation scenarios of the network, and put on a rigorous footing some previously observed numerical results.


Introduction
The collective or emergent behaviour of large networks of neurons is a topic of ongoing interest [1,2], with applications to the study of epilepsy [3,4], binocular rivalry [5], visual hallucinations [6,7] and working memory [8][9][10], among others. One type of network often considered is a ring, where neurons can be thought of as being arranged on a closed curve. This is natural if some property of a neuron is correlated with an angular variable such as heading in a head directional network [11], the direction in the plane to a visual stimulus that is to be remembered [12] or the orientation of a neuron's receptive field [13]. A variety of model neurons have been considered when studying ring networks including leaky integrateand-fire [9,10], quadratic integrate-and-fire (QIF) [14][15][16][17], more realistic conductance-based models [5,18,19] and theta neurons [20,21]. Note that under a simple transformation, the quadratic integrate-and-fire neuron with infinite threshold and reset is exactly equivalent to a theta neuron. The theta neuron is the normal form for a saddle-node-on-an-invariant-circle (SNIC) bifurcation [22,23], and we will consider networks of theta neurons  because they have remarkable mathematical properties, which allow us to perform many calculations exactly. Our results should be broadly applicable to ring networks of other Type-I neurons, i.e. neurons that start firing through a SNIC bifurcation.
Synaptic coupling between neurons is typically of longer range than just nearest-neighbour, so ring networks of neurons often have non-local coupling. This coupling is normally homogeneous, i.e. the strength of connections between two points on the ring depends on only the distance between these points [6].
There are several types of solutions of interest for a ring network of neurons. One is a 'bump' state, in which a spatially localized group of neurons is active, with all other neurons in the domain quiescent [10,24,25]. Because of the homogeneous coupling, such bumps can be positioned anywhere on the domain, and thus their position encodes a single angular variable. Such networks are often bistable, with the spatially uniform 'all-off' state being another attractor. It is this bistability that is of interest when modelling working memory [8][9][10]26].
Note that in the limit γ → 0 the distribution g(η) becomes a delta-distribution, therefore in this case we obtain system (1.1) for identical neurons. Each neuron receives a current input which is obtained as a weighted sum of pulses P n (θ ) = a n (1 − cos θ) n .
The pulse function P n (θ ) models the output pulse created when a neuron fires, i.e. θ increases through π . Note that it is possible to take the limit n → ∞, giving delta function coupling [27]. We only consider the case n = 2 below, and based on our experience and the work of others [27] we expect that varying n will only move bifurcations in parameter space, not introduce new dynamics. The parameter κ is the overall coupling strength and the weights K jk are defined according to the rule K jk = K(2π (j − k)/N) with some even non-constant 2π -periodic function K(x). Following [20], below we use a cosine coupling kernel  While neurons are either excitatory or inhibitory, connectivity functions of mixed sign are often used to mimic the combined effects of excitatory and inhibitory populations [25].
From an application point of view, for every solution of equation (1.1) it is important to clarify which neurons are quiescent and which are firing. Mathematically, this question is equivalent to calculating the firing rate for neuron k  1b. Note that the coincidence of all firing rates f k does not imply automatically the correlation of the spiking events in neighbouring neurons. In fact, the spikings can occur randomly or following a certain order such that the identical values of f k are obtained only after averaging over time. System (1.1) can also support spatially heterogeneous states, including modulated rest states, see figure 2a, and bump states, see figure 2b. Moreover, in system (1.1), one can also observe more complex non-stationary patterns, see figure 3, where more than one group of neurons may be active and their spiking events are modulated not only in space but also in time. Note that the patterns presented in figure 3 occur for parameter values for which all simpler collective states in system (1.1) are unstable. While interesting in their own right, in this paper, we study the existence and stability of only stationary patterns. The structure of the rest of the paper is as follows. In §2 we present the equations governing the continuum limit of (1.1), in §3 we discuss invariant sets of the continuum equations and in §4 we derive self-consistency equations whose solutions describe all stationary states. In §5 we show how to calculate the stability of any steady state and then in § §6 and 7, we use these results to analyse spatially uniform and non-uniform states, respectively. We conclude with a discussion in §8.

Continuum limit
Using well-established techniques, namely the Ott/Antonsen ansatz [28,29], it was shown in [20] that in the continuum limit (N → ∞), the long-term coarse-grained dynamics of system (1.1) can be described by a mean-field equation  x k x k 2 and H n (z) = a n ⎡ ⎣ C 0 + n q=1 C q (z q + z q ) ⎤ ⎦ is a nonlinear function with While the complex quantity z may not seem to have an obvious biological interpretation, one can use the equivalence of theta and QIF neurons to determine relevant quantities. Defining W ≡ (1 − z)/(1 + z) one can show [21,30] that the instantaneous firing rate of neurons at position x and time t is the flux through θ = π : If z(x, t) does not depend on t, then formula (2.2) yields a continuum limit analogue of the averaged firing rates defined by equation (1.3). Similarly, if V j = tan (θ j /2) is the voltage of the jth QIF neuron, the mean voltage at position x and time t is given by Im W.
Remark 2.1. It follows from the above definition that H n (1) = 0 for every n ∈ N. Indeed, using the substitution m = k − l, we obtain Therefore, H n (1) a n =

Remark 2.2.
Note that the function H n (z) is real even for a complex argument z. To stress this fact, we can rewrite this function in the form H n (z) = a n C 0 + 2Re D n (z) where D n (z) = a n n q=1 C q z q .
where J(x, t) = η 0 + κKH n (z(x, t)). Then, for every point x 0 ∈ [0, 2π ] it can be considered as an ODE with respect to z(x 0 , t) with a driving force J(x 0 , t). This fact and the special structure of equation (2.4) have important consequences for the dynamics of solutions to equation (2.1), which we discuss in the next section.

Invariant sets of equation (2.1)
Let us consider an ordinary differential equation with a complex-valued unknown function u(t). In the following, we suppose that J(t) is a real coefficient and γ ≥ 0. Moreover, we denote by D = {u ∈ C : |u| < 1} the open unit disc in the complex plane, by D = {u ∈ C : |u| ≤ 1} the closure of D, and by S = {u ∈ C : |u| = 1} the boundary of D. Proof. The complex conjugate of equation (3.1) reads Then, simple calculations yield For |u| = 1, this equation implies Hence, if |u(0)| ≤ 1, then |u(t)| cannot grow above one, and therefore |u(t)| ≤ 1 for all t > 0.
On the other hand, if γ > 0, then for every u ∈ S except u = −1, the differential inequality (3.2) becomes a strict inequality Moreover, at u = −1, the right-hand side of equation (3.1) equals −2i, and therefore u = −1 is not an equilibrium of equation (3.1). From all these facts, it follows that any initial condition u(0) ∈ S is pushed inside the disc D for t > 0. Moreover, if u(0) ∈ D, then u(t) also remains inside D for all t > 0. The last assertion of the proposition follows from a simple observation that for all |u| = 1 and γ = 0.

Self-consistency equation
We will start with an auxiliary proposition.

Proposition 4.1.
For any c ∈ R, the equation has exactly one stable equilibrium in the closed unit disc D. More precisely, if γ = 0, then the equilibrium is neutrally stable, but if γ > 0 then it is linearly stable.
Proof. All equilibria of equation (4.1) satisfy a quadratic equation which can be rewritten in the equivalent form where ξ denotes one of the two distinct values of the square root of c + iγ . The rational function in the left-hand side of equation (4.2) is a particular Möbius transformation, which maps the closed unit circle D on to the closed half-plane {u ∈ C : Re u > 0}. Therefore, to get a solution u ∈ D from equation (4.2), we need to choose ξ = c + iγ with Re ξ ≥ 0. In this case, equation (4.2) yields Let us linearize equation (4.1) around the above equilibrium u. In this case, the dynamics of small perturbations v obeys a linear equation Therefore, the linear stability of equilibrium u requires Using (4.2) and (4.3), we convert the expression in square brackets as follows: Therefore, the equilibrium u is (neutrally) stable if and only if Im ξ ≥ 0.
In the appendix, it is shown that for any c ∈ R and any γ ≥ 0 there exists exactly one value ξ = c + iγ such that Re ξ ≥ 0 and Im ξ ≥ 0. This ends the proof.
Let U γ (c) be a function that determines the equilibrium u ∈ D of equation (4.1) mentioned in proposition 4.1. More precisely, where the square root c + iγ is calculated according to proposition A.1 in the appendix.

Remark 4.2.
In the degenerate case γ = 0, function U γ (c) simplifies as follows: Then, if we calculate the firing rates Suppose that a(x) is an equilibrium of equation (2.1). Then a(x) is also an equilibrium of equation (2.4) with J(x, t) = w(x) := η 0 + κKH n (a). Using proposition 4.1, we find a(x) = U γ (w(x)). Therefore, the last two formulae are consistent with each other iff w(x) satisfies the integral equation The self-consistency equation (4.5) becomes particularly simple in the case of the cosine kernel K(x). Indeed, if in this case, we look for an even solution of equation (4.5), then due to remark 2.3 we have with some coefficientsŵ 0 ,ŵ 1 ∈ R. Inserting expressions (2.3) and (4.6) into equation (4.5) and equating the x-independent terms and the terms proportional to cos x separately, we obtain a nonlinear two-dimensional system which can be solved approximately with a standard Newton's method. If we find a pair (ŵ 0 ,ŵ 1 ) T ∈ R 2 that satisfies (4.7) and (4.8), then formula (4.6) yields the corresponding solution of equation (4.5), while the formula a(x) = U γ (w(x)) yields the corresponding equilibrium of equation (2.1). Note that similar self-consistency equations have been found in similar systems [31][32][33] but those authors were studying chimera states.
Having seen how to find equilibria of equation (2.1) we now move on to determine their stability.

Stability of equilibria of equation (2.1)
Suppose that z = a(x) is a stationary solution of equation (2.1). Then, inserting the ansatz z = a(x) + v(x, t) into equation (2.1) and linearizing the resulting relation with respect to small perturbations v(x, t), we obtain a differential equation determining the linear stability of a(x) and To investigate the decay of different spatial modes, we insert the ansatz into equation (5.1) and equate separately the terms at e λt and e λt . Thus, we obtain a spectral problem Note that the right-hand side of equation ( in the context of studying the chimera states [34]. Using the results of these works, we can draw the following conclusions: (i) The spectrum determined by equation (5.2) lies in a bounded region of the complex plane.
(ii) The spectrum consists of two parts: essential σ ess and discrete σ discr spectra.
(iii) The essential spectrum σ ess is determined by the multiplication operator in the right-hand side of (5.2). It can be computed explicitly by (iv) The discrete spectrum σ discr consists of a finite number of eigenvalues with finite multiplicity.
is a solution of the self-consistency equation (4.5). Then the essential spectrum corresponding to a(x) is given by Proof. Replacing a(x) by U γ (w(x)) in the definition of μ(x) and using equation (4.5), we obtain

Remark 5.2.
Recall that the branch of the complex square root in (4.4) is chosen so that Im ( w(x) + iγ ) ≥ 0. Therefore, Re (2i w(x) + iγ ) ≤ 0 and hence the essential spectrum σ ess in proposition 5.1 is always neutrally stable. Moreover, for γ > 0 the value of w(x) + iγ cannot be real, therefore in this case, the essential spectrum σ ess lies entirely in the open half-plane {λ ∈ C : Re λ < 0}.
The spectral problem (5.2), in general, is infinite-dimensional. However, in the case of cosine kernel (1.2), it can be reduced to a finite-dimensional form. Indeed, remark 2.3 implies that where V k ∈ C 2 , and φ 1 (x) = 1, φ 2 (x) = cos x, φ 3 (x) = sin x are three linearly independent functions that span the range of the operator K. Inserting (5.3) into equation (5.2) and expressing (v + , v − ) T from the resulting equation, we obtain where is a (2 × 2)-matrix. Now, using remark 2.3 and the definition of functions φ k (x) given after (5.3), we conclude Because of the linear independence of φ j (x), this formula agrees with (5.3) if and only if (Note that in the above expressions the averaging over x is performed for each matrix element separately.) Finally, we see that system ( In the next section, we consider the existence of spatially uniform steady states of equation (2.1) and use the results in this section to determine their stability.

Spatially uniform states (a) Existence of uniform states
It is easy to verify that for any cosine kernel (1.2) (in fact, even for any kernel that satisfies 2π 0 K(x) dx = 2π ) equation (4.5) has a set of spatially uniform solutions where the constant p satisfies p = η 0 + κH n (U γ (p)), or equivalently η 0 = p − κH n (U γ (p)). (6.1) By construction, formula (6.1) yields a parametric representation of all uniform states of equation (2.1) with stable essential spectra. Indeed, for every (p, κ) ∈ R 2 , we can find the corresponding excitability parameter η 0 using equation (6.1) and the corresponding equilibrium of equation (2.1) by using the formula U γ (p). Importantly, it follows from remark 4.2 that for γ = 0 the sign of the parameter p allows us to distinguish between the uniform spiking states (p > 0) and the uniform rest states (p < 0). Moreover, for small values of γ , the same criterion can be used to distinguish between the partially synchronized spiking states and the partially synchronized rest states introduced in [23]. Formula (6.1) gives a global representation of the branch of uniform states. However, in practice, it is more important to answer the other question: How many different uniform states can be found in equation (2.1) for given parameters (η 0 , κ) ∈ R 2 and what are their types? This question can be answered geometrically, if we rewrite equation (6.1) in the form where F n,γ (p) = H n (U γ (p)) = a n C 0 + 2Re (D n (U γ (p))).
Note that the right-hand side of equation (6.2) does not depend on either η 0 or κ. Moreover, for any fixed γ ≥ 0 we have U γ (p) → −1 as p → ±∞. This implies that F n,γ (p) → H n (−1) as p → ±∞, and therefore for any κ ∈ R and any η 0 ∈ R, equation (6.2) has at least one real solution. In the next proposition, we show that under certain conditions equation (6.2) can also have three different roots.  (e) there exists p min < 0 such that F n,0 (p) < 0 for p < p min and F n,0 (p) > 0 for p min < p < 0.
Proof. Suppose κ > 0. Then the existence of a unique value p 0 > 0 such that F n,0 (p 0 ) = 1/κ follows from the assumptions (a)-(c). The remaining part of the assertion (i) follows from the geometric consideration of equation (6.2) explained in figure 4. Indeed, if the slope κ is fixed, then changing η 0 we can encounter three qualitatively different situations, which are separated by two critical cases: the line (p − η 0 )/κ is tangent to the graph of F n,0 (p) for p > 0, or this line passes through the point (0, 0).

Remark 6.3.
It is easy to verify that the assumptions of proposition 6.1 are satisfied for n = 2. In particular, The conclusion of proposition 6.1, in this case, is represented geometrically in figure 5.

(b) Stability of uniform states
Let w(x) = p 0 ∈ R be a constant solution of equation (4.5). It corresponds to the equilibrium z = U γ (p 0 ) of equation (2.1), and therefore we can perform its stability analysis according to the approach described in §5. Using proposition 5.1, we find that the essential spectrum σ ess consists of two points μ 0 and μ 0 where μ 0 = 2i p 0 + iγ . To calculate the discrete spectrum σ discr , we need to consider equation (5.6). Note that the matrix L(x, λ) of a uniform state does not depend on x. Moreover, simple calculations yield Note that equations (6.5) and (6.6) obviously coincide with each other and therefore determine identical eigenvalues.

Proposition 6.4. Equation (6.4) has two solutions
while equation (6.5) has two solutions Proof. Let us consider equation (6.5). Using the definition of B jk (λ), we rewrite it in the form For every λ / ∈ σ ess (that means λ = μ 0 and λ = μ 0 ), the latter equation is equivalent to Solving it with respect to λ we obtain the formula for λ 2,± . The other equation (6.4) can be considered by analogy.
Remark 6.5. Using (4.4), it is easy to demonstrate that that determines λ 1,± and check when it has at least one zero solution. This situation corresponds to a static bifurcation of the uniform state a 0 = U γ (p 0 ) and is equivalent to the condition Inserting here the expression of ζ 0 from remark 6.5, we obtain The case μ 0 = 0 corresponds to the degenerate essential spectrum and therefore we discard it. Then the only possibility to satisfy the above equation is to have It is easy to verify that the latter condition is equivalent to the situation when the line in the lefthand side of equation (6.2) is tangent to the graph of F n,γ (p). In other words, all static bifurcations described by equation (6.4) occur at the fold points of the graph of p versus η 0 with fixed κ. Proposition 6.6. In the case γ = 0, the eigenvalues λ 1,± determined by equation (6.4) are either both real or complex-conjugate and purely imaginary. Therefore equation (6.4) in this case can determine only static bifurcations, i.e. bifurcations characterized by real positive eigenvalues appearing from zero.
Proof. If p 0 > 0, then μ 0 = 2i √ p 0 and the derivatives U 0 (p 0 ) and D n (U 0 (p 0 )) are both real. Hence, the value of ζ 0 is purely imaginary and therefore On the other hand, if p 0 < 0, then μ 0 = −2 √ −p 0 and therefore This ends the proof.
See figure 6 for a graphical representation of this proposition.
Proof. The existence of the uniform state z = U γ (p 0 ) is obvious and we need to consider only the corresponding expression taking into account that μ 0 and ζ 0 depend on γ . Due to the assumptions p 0 > 0 and 1 − κF n,0 (p 0 ) > 0 we can be sure that the square root in equation (6.7) remains purely imaginary for sufficiently small γ , therefore Re λ 1,± = Re (μ 0 + 2ζ 0 ).

Spatially non-uniform states
In the previous section, we considered the constant equilibria of equation (2.1), which correspond to uniform rest states and uniform spiking states of coupled theta neurons (1.1). However, these are in general not the only possible solutions of equation (2.1). Below we show that for a cosine coupling (1.2) with A = 0 equation (2.1) also has a variety of non-constant equilibria representing more complex patterns in system (1.1). We focus on the case γ = 0, n = 2 and test three distinct values κ = 1, −1 and −2 corresponding to three qualitatively different scenarios in figure 5. For each of these cases, we consider three distinct values of A, namely 0, 3 and −5.
In figure 7, for all nine combinations of κ and A values, we plot either p (for spatially uniform states) orŵ 0 (for non-uniform states) as a function of η 0 . The branch of uniform states is always independent of A, but the stability of different states on it can change for different A. If A = 0, then using proposition 6.4, we typically find two bifurcation points determined by the condition λ 2,± = 0. Then, using the self-consistency equations (4.7)-(4.8) we can compute a branch of nonconstant equilibria of equation (2.1) that connects these two bifurcation points. Every solution (ŵ 0 ,ŵ 1 ) of equations (4.7)-(4.8) withŵ 1 = 0 corresponds to a non-constant even solution (4.6) of equation (4.5) and hence to a spatially modulated equilibrium a(x) = U γ (w(x)) of equation  2p 2p (2.1). The stability of this equilibrium can be analysed using the results of §5. More precisely, the essential spectrum σ ess of the operator on the right-hand side of the linearized equation (5.1) is determined by proposition 5.1, while the discrete spectrum σ discr of this operator can be calculated by numerically solving the characteristic equation (5.6). Note that since w(x) is an even function, the matrix L(x, λ) in the definition of B(λ) is an even function of x too. Therefore, B jk (λ) = 0 for j = 1, 2 and k = 3 as well as for j = 3 and k = 1, 2. Consequently, equation (    Each row of these figures shows an equilibrium a(x) of equation (2.1), the corresponding firing rate f (x) calculated by (2.2), and the linear stability spectrum of a(x). The sampled equilibria are chosen so that they illustrate all qualitatively different spatial profiles a(x) and their linear stability properties. Note that we do not show any equilibria for figure 7i, because they look very similar to those shown in figure 12.
Our results reveal several parameter combinations suitable for the observation of stable bump states. In particular, such states can be found for (κ, A) = (1, 3), (κ, A) = (−1, −5) and (κ, A) = (−2, −5). Although in all these cases, the bump states have similar spatial profiles, their stability diagrams for positive and negative κ are a bit different. More precisely, for κ = 1, the η 0 -interval of stable bumps is bounded by two fold points, while for κ = −1 (as well as for κ = −2) the stability interval of bumps is bounded by a fold point at the right end and by a Hopf bifurcation at the left end.
Apart from the bump states, equation (2.1) has only one more type of stable spatially nonuniform equilibria. These are modulated rest states shown in the first and third rows of both figures 10 and 11. Even though none of the neurons are firing here, their state is not spatially uniform.
Finally, we emphasize that our consideration reveals also a wide parameter range, see figure 7c,e,f,i, where all equilibria of equation (2.1) are unstable. In this case, any simple collective dynamics of neurons is not possible, therefore more complex non-stationary patterns can emerge in system (1.1), see figure 3.

Discussion
In summary, we analysed the dynamics of a ring network of theta neurons, homogeneously coupled with a cosine kernel. We took the continuum limit that allows us to describe the network's asymptotic dynamics by a complex-valued integrodifferential equation (2.1). We showed how to determine the existence and stability of stationary solutions, both spatially uniform and nonuniform. We have given a number of examples demonstrating the validity of our approach, mostly concentrating on the case of identical neurons. These results put on a rigorous basis some of the numerical observations in [20,21]. We now discuss how our results relate to those of others. Byrne et al. [17] considered either one or two populations of theta neurons with synaptic coupling using reversal potentials and characteristic timescales, and exponentially decaying coupling kernels. They considered domains with either periodic or Neumann boundary conditions and looked for Turing bifurcations of the spatially uniform state. They also performed numerical bifurcation analysis of time-dependent solutions such as wavetrains and travelling fronts. Esnaola-Acebes et al. [14] also considered two ring networks of excitatory and inhibitory QIF neurons with even coupling functions involving more than just the first cosine term. They chose parameters so that there was a single spatially uniform state and performed a Turing bifurcation analysis of this state. This analysis was used to explain the decaying oscillatory transient spatial modes seen when the spatially uniform state was perturbed. They also numerically followed a bump solution and found a scenario of the form shown in figure 7e. Byrne et al. [16] considered a network of QIF neurons on a periodic domain, with synaptic dynamics and an inverted Mexican hat connectivity. They also included gap junctional coupling. By analysing the continuum equations they found Turing-Hopf bifurcations of the spatially uniform state that result in the appearance of travelling and standing waves. Schmidt & Avitabile [15] studied ring networks of QIF neurons with Mexican hat connectivity and also determined the conditions for a Turing bifurcation.
Most of these works presented largely numerical studies of the fully nonlinear and often timedependent solutions that such networks can support whereas we have concentrated on stationary solutions, exploiting the form of the coupling function to explicitly calculate these and analytically determine their stability.
Our approach could easily be generalized to other forms of coupling such as synaptic coupling using reversal potentials [17] or gap junctional coupling [16,35], and different coupling kernels, as long as the form of the equations allows the use of the Ott/Antonsen ansatz.
Data accessibility. The codes used for producing the data for figures are available in electronic supplementary material [36].